﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "pqinteg")]
    public static unsafe double pqinteg(double a0, double b0, double eps, IntPtr f_x_ptr)
    {
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return pqinteg(a0, b0, eps);
    }

    /// <summary>
    /// 连分式求积法
    /// f计算被积函数值f(x)的函数名。
    /// </summary>
    /// <param name="a0">积分下限。</param>
    /// <param name="b0">积分上限。要求b0>a0。</param>
    /// <param name="eps">积分精度要求。</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double pqinteg(double a0, double b0, double eps = 1.0E-7)
    {
        int k, j, il, flag, n;
        double h0, g0, d, s0, s1 = 0.0, x;
        double* b = stackalloc double[10];
        double* h = stackalloc double[10];
        double* g = stackalloc double[10];

        il = 0;
        n = 1;
        h0 = (b0 - a0) / n;
        flag = 0;
        // 梯形公式计算初值
        g0 = h0 * (f_x(a0) + f_x(b0)) / 2.0;
        while ((il < 20) && (flag == 0))
        {
            il = il + 1;
            h[0] = h0;
            g[0] = g0;
            // 计算b[0]
            b[0] = g[0];
            j = 1;
            s1 = g[0];
            while (j <= 7)
            {
                d = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    x = a0 + (k + 0.5) * h[j - 1];
                    d = d + f_x(x);
                }
                // 变步长梯形求积法计算新近似值g[j]
                g[j] = (g[j - 1] + h[j - 1] * d) / 2.0;
                h[j] = h[j - 1] / 2.0;
                n = 2 * n;
                // 计算b[j]
                funpqj(h, g, b, j);
                s0 = s1;
                // 连分式法计算积分近似值s1
                s1 = funpqv(h, b, j, 0.0);
                if (Math.Abs(s1 - s0) >= eps)
                {
                    j = j + 1;
                }
                else
                {
                    j = 10;
                }
            }
            h0 = h[j - 1];
            g0 = g[j - 1];
            if (j == 10)
            {
                flag = 1;
            }
        }
        return (s1);
    }

    /*
    // 连分式求积法例
      int main()
      { 
          double d,eps,pqintegf(double);
          eps=0.0000001;
          d=pqinteg(0.0,4.3,eps,pqintegf);
          cout <<"积分值 s = " <<d <<endl;
          return 0;
      }
    // 计算被积函数值
      double pqintegf(double x)
      { 
          return(exp(-x*x));
      }
    */
}

